Case study of Rb + (aq), quasi-chemical theory of ion hydration, and 
the no split occupancies rule 
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Quasi-chemical theory applied to ion hydration combines statistical mechanical theory, electronic structure calculations, and 
molecular simulation, disciplines which are individually subjects for specialized professional attention. Because it combines 
activities which are themselves non-trivial, quasi-chemical theory is typically viewed with surprise. Nevertheless, it provides a 
fully-considered framework for analysis of ion hydration. Furthermore, the initial calculations are indeed simple, successful, and 
C") provide new information to long-standing experimental activities such as neutron diffraction by hydrated ions. Here we review 
,— ' quasi-chemical theory in the context of a challenging application, Rb + (aq). 
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Water is a chemically active liquid. Probably the most prim- 
itive aspect of that chemical activity is dissolution of elec- 
trolytes, and the chemical processes based on availability of 
dissolved ions. One example is salt as a trigger of autoim- 
mune diseased™ Another example, presumably related to the 
frrstj^is the selective transport of dissolved ions across mem- 
branes. 7 " 9 Rubidium (Rb + ) is interesting in this respect be- 
cause it serves as an analog of potassium (K + ) that conducts 
current through potassium ion channels, even though Rb + is 
slightly larger (by 0.2 AIDSES 

In addition to this chemical activity, water has long been a 
serious challenge for statistical mechanical theory of liquids, 
which itself is properly almost entirely classical mechanical 
theory. 12 The challenge presented by liquid water is the vari- 
ety of intermolecular interactions that must be considered with 
a wide range of interaction strengths (Fig.JTJ.ElThose interac- 
tions include excluded-volume repulsions, essential since liq- 
uids and liquid water are dense materials. Those interactions 
also include H-bonding interactions that are attractive on bal- 
ance, much stronger than thermal energies, and essential for 
the characteristic behavior of liquid water. 

Attractive interactions that involve many neighbors are 
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Fig. 1 Probability densities for pair contributions to the binding 
energy of a water molecule to liquid water.E3 Rightmost curve: the 
most positive (unfavorable) pair contribution; leftmost curve: the 
most negative (favorable) contribution. The middle curve (with 
yellow-shaded area) is the probability density of all the pair 
interaction contributions.^ 



good cand idates for treatment by mean-field approxima- 
tions.^^ Of course, if those attractive interactions are indi- 
vidually weak on a thermal energy scale, that characteristic is 
favorable for simple theories also. 

Here we consider rubidium dissolved in aqueous solution, 
Rb + (aq), to provide a contrast with theories of liquids more 
broadlyQ^GIIand to pursue a specific discussion of what is yet 
required theoretically to treat solvated electrolytes. The hydra- 
tion free energy of Rb + (aq) is known to be roughly —100^8^ 
under standard conditions, favorable enough to dissolve sim- 
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pie Rb + salts, and indeed large on the thermal energy scale. 
Additionally (Fig.|2|, the number of near-neighbors is modest, 
between four and seven. As another contrast, t he iso electronic 
Kr(aq) has about eighteen (18) near-neighbors.l^^From this 
comparison it is clear that the features of strong attractive in- 
teractions and a reduced number of near-neighbors are corre- 
lated: the reduced number is a consequence of the crowding 
of near-neighbors drawn close by the interesting attractive in- 
teractions. 

This discussion suggests that we seek a way forward by 
focusing on the small number of near-neighbors, deploying 
direct quantum mechanical computation for Rb(H20)„ + with 
a small number n of neighbors, then stitching those compu- 
tational results into the broader theory of liquids. That was 
indeed the idea of quasi-chemical theory discussed here. At 
a formal level that theory is fully conclusive, but it is sur- 
prising that it was not worked-out until fairly recentlyJ^QHSl 
It was also surprising that the initial applications of the theory, 
to Li + (aq), 26 were highly effective - predictions of hydra- 
tion free energy matched experiment 27 and ab initio molec- 
ular simulation estimates. 28 Further, the initial applications 
provided new information to long-standing neutron diffraction 
work on hydrated ions, 29 and indeed motivated renewed effort 
on those experiments.^ 

The initial applications of quasi-chemical theory could be 
understood on a physical and intuitive basis. That encouraged 
the generation of simple mimics that were less fully thought 
through. A focused discussion of the statistical thermody- 
namic formalities was given by Asthagiri, et ai, 31 with the in- 
tention of reducing the confusion that can result from a crowd 
of imitators. The presentation of Asthagiri, et a/.^Hwill be the 
basis of the discussion that follows. 

2 Theory 

The primary theoretical target for quasi-chemical theory has 
been the excess (or interaction part) of the partial molar Gibbs 
free energy (or chemical potential) of the species of interest, 
here the Rb + (aq). This is indeed a basic characteristic of 
the solution and the Rb + ion in it, but it is also comparatively 
simple. The potential distribution theorem 12 offers a partition 
function for evaluation utilizing information obtained on the 
local environment of the ion. It is found that the desired free 
energy can be cast as 




Fig. 2 Ab initio molecular dynamics (AIMD) result for the radial 
distribution of water oxygens about Rb + , gRbo- We used 
VASP-version 4.2pEIwith the Perdew- Wang (PW91) 
exchange-correlation functional,! 34 ! 35 ! a plane-wave basis set cutoff 
at 36.75 Ry with the interaction between valence and core electrons 
described by the projector augmented-wave method (PAW).I 3 -^ 
Corrections to capture long-ranged interactions like dispersion have 
not been included, and recent studies of ions in water suggest that 
they are not helpful. 37 The system consisted of one Rb + ion and 64 
waters in a cubic simulation cell with edge length of 12.417 A and 
full periodicity. Charge balance was achieved with a neutralizing 
background. All hydrogen atoms in the system were replaced by 
deuterium. The simulation was carried out for 49.28 ps with a time 
step of 0.5 fs. The average temperature during the simulation was 
347.9 ±4.5 K, set intentionally high to avoid the over- structuring of 
pure water observed in room temperature AIMD simulations J 38 " 41 1 
The numbered curves show a neighborship analysis from 
decomposing £ri,o - curves in red highlight waters that occupy an 
intermediate region distinct from the first «=l-4 waters that fill in 
the peak of the first maximum, and more distant waters that mostly 
occupy the second maximum. Results over longer length-scales are 
shown in Fig. [4] . 



treatment of the complex itself as a chemical cons tituent is a 
characteristic feature of this quasi-chemical theory El 22 l 24 l 25 1 
Other features of Eq. ([T]) properly fill-in a picture of the 
chemical equilibrium of this complex. The combination 
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sociation reaction 



n[i H q begins a free energy balance for the as- 



Rb+ +«H 2 # Rb(H 2 0)„~ 



(2) 



A t Rb+ = -kTlnK^p H2 o n +kTlnp Rb+ («) 

,(ex) 



^Rb(H 2 0),t-^H 2 0- (1) 



On the right side, notice the reference to the molecular com- 
plex Rb(H20) n + and its excess free energy M R b(H,o) +' me 



Similarly, 



K„ 



(0) _ PRb(H 2 0)„-t 



P Rb +PH 2 0" 



(3) 



with px the number density of species X, is the equilibrium ra- 
tio for that association reaction (Eq. (|2|) treated as in an ideal 
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Fig. 3 The coordination number distribution, for the inner-shell 
radius of A = 3.76 A, from the AIMD simulations and an "LJ 
model." For the LJ model, we used Aqvist Lennard- Jones 
parametersEHfor the Rb + ion, (e = 1.71 x 10~ 4 kcal/mol, 
(7 = 5.62177 A), and the SPC/E potential 43 for water intermolecular 
interactions. We carried out standard NVT molecular dynamics 
calculations using the GROMACS 44 package (version 3.1.4). The 
system consisted of one Rb + and 2177 water molecules in a (40 A) 3 
cell. A single CP was included for charge balance. For electrostatic 
interactions, the particle mesh Ewald technique was implemented 
with Fourier spacing of 1.5 A, a sixth-order interpolation, a 10 A 
cutoff in a direct space, and a tolerance of 10~ 5 . A cutoff distance of 
16 A was adopted for the Lennard-Jones interactions. 
Intramolecular geometric constraints on water molecules were 
enforced by the SETTLE algorithm. 45 Data was collected during a 
1.0 ns production phase, following a 1.0 ns equilibration phase. The 
Nose-Hoover thermostat with a coupling constant of 0.2 ps 
maintained T = 298.15 K J 46 " 48 1 The variation in coordination 
number n is substantially smaller for the LJ model data than for the 
AIMD results; but remember that the AIMD results correspond to 
the slightly higher T & 350 K. 



gas, i.e., neglecting interactions with the solution external to 
the complex; that is signified by the superscript zero. As em- 
phasized previously,! 12 l 22 l 3 1 1 K„ is a well-defined few-body 
computational target, and the full force of available quan- 
tum mechanical computational methods can be brought to 
bear. This has been especially relevant to treatment of transi- 
tion metal ions where d-orbital splitting may be addressed,^ 
and where the simplest applications of the present theory are 
particularly effective. 49 50 Electronic charge-transfer between 
ion and ligands is another effect of quantum-mechanical ori- 
gin that is transparently and simply included in these quasi- 
chemical approaches. Proper accounting of electronic charge 
density is important for computing absolute and relative ion 
binding affinities, which is useful for probing mechanisms of 
selective ion binding.EO 



The remaining feature of Eq. ([T]) is the probability p Rb + (n) 
of observing n ligands within a defined inner shell (Fig. [3]). 
That inner shell is a fundamental concept for this approach, 
and we return below to discuss it further. For now, note that 
if only one coordination number n were ever observed, then 
p Rb + (n) = 1 and that contribution in Eq. ([1} would vanish. 
Ultimately, that contribution in Eq. ([TJ carries the full ther- 
modynamic effect of whatever actual variablity of inner-shell 
occupancy does occur. 

Leaving the occupancy probability for general consider- 
ation, the left side of Eq. ([TJ is independent of n, and 
Eq. ([T| therefore describes the n-dependence of the probability 
Pm>+ ( n )- If our 8 oa l i s to evaluate the free energy, however, 
we can choose n for our convenience. An interesting choice is 
n = h, the most probable value (Fig.|4jl. This choice makes the 
negative contribution kT In p Rb + (n) as small as possible, and 
suggests neglecting that contribution to obtain the convenient 
approximation 

The neglected contribution is negative, and the approximate 
result Eq. Q will be higher than the true free energy. Nev- 
ertheless, a specific evaluation of p Rh + (h) can be straight- 
forwardly extracted from standard molecular simulations,! 52 ! 53 ! 
and thus the significance of fluctuatations of composition of 
the inner shell is obtained merely by noting the size of the 
neglected kT In p Rb + (it). Those results are obtained and dis- 
cussed below. 

3 Application to Rb + (aq) 

The primitive quasi-chemical theory, Eq. (p), has indeed been 
applied to several hydrated metal ions 2 "! 50 ! 54 " 57 ! as well as 
other solvation problems.^HH For transition metals, as ex- 
amples, the near-neighbor water molecules are clearly located 
on the basis of chemical considerations, 49 and this theory is 
straightforwardly successful. Generic procedures for those 
standard cases were given by Pratt and Asthagiri. 24 For other 
cases, some physical judgement is required, and evaluation of 
the various contributions to Eq. ([1} requires analysis. We con- 
sider application to Rb + (aq) to show how that goes. 

3.1 No split occupancies 

It is important that theories teach how to understand physi- 
cal problems in addition to reproducing numerical values of 
central properties. For our present problems, that learning is 
focused on how to catagorize near-neighbor water molecules 
to make simple theories effective. 

For cases like Rb + , in contrast to transition metals, the im- 
portant observations arise from the AIMD simulation results 
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Fig. 4 Distribution gRboM of oxygen (O) atoms radially from the 
Rb + ion in the AIMD simulation, emphasizing the variation of the 
most probable occupancy (n, the blue dots) of the spherical inner 
shell defined by the indicated radius r. n, being integer-valued, 
changes discontinuously with defining radius, but reliably tracks the 
mean number of oxygen atoms within r, i.e., (n(r)} = 
4^po Jo gRboM^d* f° r water at density Pq. For example, (n(r)) = 
6.86 ± 0.21 and h = 7 at r = 3.76 A for AIMD compared to 7.43 ± 
0.22 and 7 from the LJ model. Nevertheless, note that (n(r)) does 
not exhibit a convincing plateau for physical identification of a 
coordination number because the principal minimum of gRbo near r 
= 3.76 A is weak. 



of Fig. [2] The first minimum of the radial distribution is re- 
markably mild and does not provide a convincing identifica- 
tion of an inner shell. Instead, we consider the neighborship 
decomposition of that radial distribution. We see that the 7th- 
most nearest neighbor (n=7) contributes to the first maximum 
(with negligible contribution to the peak), the second maxi- 
mum, and the first minimum. In other words, the contribution 
of the 7th-most distant water neighbor to g(r) is multi-modal. 

To make our problem simple, we try to set an inner-shell 
volume so that the 7th neighbor does not contribute. We see 
from Fig.|5]that p Rb + (n = 7) is particularly small for an inner- 
shell radius of X = 3.2 A. With the indicated choice of inner 
hydration shell, we notice further that neighbors 1-4 fill-out 
the principal maximum of that radial distribution function. 

This lesson we will call the no split occupancies rule. 
These neighborship analyses^have become characteristic of 
quasi-chemical theor ies and were used previously for Li + (aq), 
Na+(aq) and K+Caq).' 25 ! 29 ' 56 ' 5 ^ 



3.2 Methods 

To evaluate the quasi-chemical free energy contributions to 

(ex) \ 
^Rb+ ' We Start W ^ t ^ le ^ rSt term °^ ^ ^ n ' S term §i VeS 

the free energies for association of Rb + with n water ligands 
to form clusters within our choice of inner-shell radius (X = 
3.2 A). These clustering equilibria take place in an ideal gas 
with a water density corresponding to a pressure of 1 atm. 

Gas-phase thermochemical data required for the association 
reactions (Eq. |2|) were obtained by electronic structure cal- 
culations using the Gaussian09 program 64 and density func- 
tional theory with Becke's three-parameter exchange func- 
tionallSland the LYP 66 electron correlation function (B3LYP). 
All structures were fully optimized with a basis including po- 
larization and diffuse functions (6-31 1++G(2d,p)) on oxygen 
and hydrogen centers, and the LANL2dz effective core po- 
tential and basis set on Rb + . At the lowest-energy geometry, 
confirmed by zero-valued slopes of the electronic energy with 
respect to atomic displacements, a standard Hessian analysis 
was performed to compute normal mode vibrational frequen- 



cies 1 1 'using the same basis set. Quantum mechanical parti- 
tion functions were then calculated, 69 thus providing a deter- 
mination of the free energy changes of the association reac- 
tions due to atomic motions internal to water and the clusters 
at temperature T = 298 K and 1 atm pressure. 

In a subsequent step, the cluster results were adjusted with 
a ligand replacement contribution, ralnpH 2 o> to account for 
the actual concentration of water ligands in liquid water at the 
density Ph 2 o =1 g/cm 3 . If this density is tracked as an adjust- 
ment of the ideal gas pressure, then it corresponds to a pressure 
factor of 1354 atm. 



To compute the last two terms of Eq. (HI), jj. 



(ex) 

Rb(H 2 0) + 
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Fig. 5 Coordination number distributions within inner hydration 
shells defined by various X: upper AIMD, lower LJ model. The 
AIMD results display a modest but distinct tendency toward lower 
coordination numbers. This is possibly connected to the generally 
incorrect description of overlap repulsions by LJ interactions. 



n/ijj q, we treated the solvent external to individual water lig- 
ands and the inner-shell clusters as a reaction field using a po- 
larizable continuum model (PCM) J- 7 -^ We subtracted the gas- 
phase electronic structure energy for the n-coordinate cluster 
geometry to obtain the desired excess free energies. These two 
terms combined make up the outer-shell electrostatics con- 
tribution. Finally, we evaluated Rb + hydration free energy, 

(ex) 

/i^ b +, by summing the quasi-chemical components for forma- 
tion of the most probable complex (Eq.|4ji, Rb(H20)^ + . 
The integral equation formalism (IEF-PCM) was imple- 

I79I 

mented for the outer-shell elecrostatics calculations. 1 ^ A ra- 
dius of 3.2 A around Rb + defined the inner-shell boundary. 
In addition, default parameters were used to define hydrogen 
and oxygen radii used to create the solute cavity as a set of 
overlapping spheres. The dielectric constant of the outer-shell 
medium was set to 78.35 to represent liquid water. 

4 Results 

The principal maximum of the Rb-O radial distribution func- 
tion (Fig. |2]i is near r ss 3.0 A. Experimental studies re- 
port similar results for the location of this maximum:l21E3 
2.93 ±0.3 A, 2.90 A, and 3.05 A. 

The free energy results (Fig. [6]i show that the coordina- 
tion number n = 4 is indeed the most probable within the in- 
ner hydration shell defined by A = 3.2 A. This then implies 
( a l) = — 65.4kcal/mol, which agrees reasonably with the 
experimental value of —69.51 kcal/molJ^S 

Addressing Eq. Q, several approximations have accumu- 
lated. The first of those is the neglect of population fluctuation. 
From Fig. |5] we see that this error amounts to kT\n p Rb + (re) w 
—0.6 kcal/mol, roughly a 1% error on the predicted hydration 
free energy. This could be easily appended to the final result, 
but it is not significant here. 

Further approximations entered to evaluate the free ener- 
gies on the right side of Eq. Q. A normal mode analyis 
yields harmonic frequencies that are expected to represent the 
vibrational motions for small ion-water clusters. A perturba- 
tive analysis 76 of anharmonicity in the electronic energy sur- 
face confirmed that vibrations in the n = 4 cluster are well- 
described by normal mode analysis. If anharmonicity were 
important, the corrections could easily be included in the final 
result. 

A serious approximation is the treatment of the exter- 
nal environment as a diectric continuum when considering 
^RbfH o) + — n l J Ti*b- This approximation is clearly not real- 
istic on a molecular scale, i.e., the solvent is not actually a 
dielectric continuum. But in this application, the dielectric 
continuum model is used for outer-shell electrostatic effects, 
and thus molecular-scale inaccuracies should be less serious. 
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Fig. 6 Contributions from the primitive quasi-chemical evaluation 

(ex) 

of /x^ b+ (aq). The predicted hydration free energy is —65.4 kcal/mol 
when the outer-shell electrostatics contribution is based on a single 
n = n = 4 solvent-adapated cluster configuration. Using instead 
many configurations sampled from the AIMD simulation record, the 
outer-shell electrostatics contributions are similar, but slightly lower, 

(ex) 

and predict flLj (aq)= —70.5 kcal/mol. The modeled inner-shell 
occupancy distributions PRb+ («) resemble the observations from the 
AIMD simulation (Fig. pi X = 3.2 A), plotted with open boxes 
where the constant 'c' is adjusted to match the single-point QCT 
model at n = n = 4. This qualitatively satisfactory comparison 
achieved with both single and multiple solvent-adapted structures 
(highlighted by shading) stands in contrast with a previous primitive 
assessment)^ for K + (aq) that is expected to be physically similar. 
That previous work did not attempt a detailed examination of the 
theory, nor did it consider specificall y a neighbor ship analysis 
(Fig.pl, as has been long customary .E322EH57159] 



Moveover, because of the balance of free energies in this con- 
tribution, there is an opportunity for molecular-scale inaccu- 
racies to cancel to some extent. 

It is also remarkable and approximate that we calculated 

(ex) 

Mdu/o n \+ on tne basis of a single cluster configuration and 

Rb(H 2 U) n 

electronic charge distribution. Moreover, we based the es- 

(ex) 

timate of ul.L _. + on the geometry and electronic charge 

Kb(ti2U) n 

distribution of the cluster that has been subtly altered by the 
environment. This inclusion of a solvent reaction field is ex- 
pected 25 to be an improvement for simple dielectric and Gaus- 
sian distribution theories of solvation free energies. Further, 
this approach represents a distinct change from typical proce- 
dures used previously. 25 Although the cluster structures (n=\- 
8) are only subtly different on the basis of casual inspection, 
the results of Fig. [6] for n > 5 are decisively improved. With 
previous procedures, the present results become qualitatively 
like the published, and puzzling, results for K + with n > 5!^ 
Of course, the similarity of Rb + and K + is the natural physi- 
cal expectation. The implied computed hydration free energy 
based on n = 4 is not much changed in quality compared to 
estimates based on the earlier procedures, 25 as also noted pre- 
viously.E^For example, is only 2 kcal/mol more positive 
using Eq. |4]based on a single cluster geometry and electronic 
charge distribution independent of the environment. The dif- 
ference here is that now the overall occupancy distribution 
(Fig. |6| is also qualitatively reasonable. 

This single-point estimate of the cluster free energy was 
tested by sampling inner-shell structures from the AIMD sim- 
ulation record and evaluating electrostatic contributions to the 
free energies using a simple dielectric continuum solvation 
model. In that case, the excess free energies of the sampled 
clusters and individual water ligands were estimated using the 
APBS-version 1.3 software^ with the same parameters de- 
scribed earlier.^ Results for the sampled clusters were then 
combined in the thermodynamically consistent fashion, i.e., 
by adding the inverses of the Boltzmann factors. That test pro- 
duced slightly lower outer-shell electrostatic contributions and 

(ex) 

/i^ b +(aq)= —70.5 kcal/mol, thus confirming the results above. 
The conclusion appears to be that isolated n > 5 clusters are 
sufficiently unusual in structure to cause trouble for single- 

(ex) 

point estimates of ir based on gas-phase structures. 

Single solvent-adapted structures, or structures sampled from 
liquid-phase simulations, result in a satisfactory improvement. 



5 Concluding Discussion 

Using AIMD simulation and primitive quasi-chemical theory, 
we find that n = 4 waters preferentially solvate Rb + within a 
spherical inner shell defined by radius X = 3.2 A. This bound- 
ary extends slightly beyond the first peak of the Rb-O radial 
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distribution function observed in AIMD simulations and re- 
ported in experiments (r = 3.0 A), but lies well within the pu- 
tative minimum at r = 3.8 A occupied by n = 7 waters. As cus- 
tomary, we intentionally chose a smaller inner hydration shell 
for free energy analysis to avoid split occupancy between first 
and second hydration shells observed for the 7th-most distant 
water from Rb + in AIMD simulations. 

The hydration free energy for Rb + predicted by primitive 
quasi-chemical evaluation agrees reasonably well with exper- 
iment. Good agreement is achieved when treating the most 
probable aqueous coordination complex within the defined in- 
ner hydration shell as a single energy-optimized gas-phase 
species using quantum mechanical methods and coupling that 
cluster to an implicit model of the solution. For example, cou- 
pling Rb + (H20)4 to a dielectric continuum model, whether or 
not the cluster adapts to that environment, results in a hydra- 
tion free energy of s» —65 kcal/mol. In comparison, coupling 
many n = 4 configurations sampled from AIMD simulation to 
an implicit solvation model results in a similar hydration free 
energy prediction of -70.5 kcal/mol, a result within 1 kcal/mol 
of experiment. 

Although the overall hydration free energy remains simi- 
lar, a significant improvement does occur in the distribution 
of inner-shell coordinations predicted when treating solvent- 
adapted clusters (one or many) compared to prior procedures 
that treated only gas-phase clusters. Comparison of the oc- 
cupancy distribution from AIMD with those from a LJ model 
indicates that AIMD results exhibit slightly lower coordina- 
tion numbers with greater variability. 

Based on similarity in size (within 0.2 A) and identical 
charge, Rb + and K + ions are expected to share similar sol- 
vation characteristics in water and other environments. Fore- 
most, hydration free energies are nearly identical according to 
experiments, a result also predicted by quasi-chemical anal- 
yses. Further similarities can be highlighted by comparing 
AIMD simulation and quasi-ch emical free e nergy analysis of 
Rb + with earlier studies of K + l 2 5 | 29 | 56 | 57 | 59 | p or exam ple, me 
first peaks in ion-0 radial distibution functions fall in simi- 
lar locations, but slightly closer for the smaller K + ion (r = 
2.8 A), as anticipated. In both cases, a weak minimum ob- 
scures identification of a coordination shell. For K + , this min- 
imum ostensibly occurs closer (r = 3.5 A) with occupation by 
n = 6 waters, one less than Rb + . Similar to Rb + , split occu- 
pancy observed in the AIMD record of K + (aq) simulations, 
specifically relating to the 6th-most distant water, motivated 
definition of a more restricted inner hydration shell for free 
energy analysis. 

Perhaps unexpectedly, these inner shells are defined at 
nearly identical distances for Rb + and K + ions (X K + = 3.1 A) 
and have similar properties: the first n=l-4 waters fill-in the 
first peaks in the ion-oxygen radial distribution functions, and 
n=4 waters preferentially solvate both K + and Rb + according 



to AIMD simulation and quasi-chemical analysis. In contrast, 
n=4 waters fill-in the first peak and preferentiall y solv ate Na + 
within a closer inner- shell distance of A = 2.6 A.EE3 Finally, 
stable inner-shell hydration structures for both K + and Rb + 
exist in gas phase that are absent in AIMD simulations of the 
liquid phase. For example, the n = 8 non-split occupancy is 
a rare composition in AIMD simulation, but forms a stable 
skewed cubic structure with 4 waters in a plane above and 4 
below K+EHandRb+EH! in the absence of stabilizing interac- 
tions with the more distant solvation environment. This high 
coordination contrasts with smaller inner-shell coordinations 
(« < 6) reported for Na + PSlFurther, these 8-coordinate clus- 
ters resemble the crystallographic ligands resolved around K + 
and Rb + ions in the binding sites of potassium ion cha nnels 
that conduct K + and Rb + , but reject smaller Na + ionsJ 10 ! 7 ^ 
Fol lowing p revious studies of selective K + and Na + bind- 
j n g |9157l59160j q Uas i- C hemical theory may be useful in future 
work to analyze the subtleties of Rb + binding and conduction 
in potassium channels. 

Quasi-chemical theory applied to ion hydration combines 
statistical mechanical theory, electronic structure calculations, 
and molecular simulation, disciplines which are individually 
subjects for specialized professional attention. Because it 
combines activities which are themselves non-trivial, quasi- 
chemical theory is typically viewed with surprise. Neverthe- 
less, it provides a fully-considered framework for analysis of 
ion hydration. 

It is striking that the three sub-disciplines noted (statisti- 
cal mechanical theory, electronic structure calculations, and 
molecular simulation) are so distinct. Typical practice in each 
subdiscipline is to parameterize the ingredients from the other 
two in order to eliminate those complexities. Thus, for ex- 
ample, sophisticated electronic structure calculations are done 
with solution models (dielectric continuum models) that are 
not justified on the basis of more basic observation. Simi- 
larly, sophisticated statistical mechanical theory is typically 
pursued where molecular-scale realism of the model can be 
empirically eliminated, for example by treating parameter- 
ized pair-decomposable models of intermolecular interactions. 
Simulation calculations also, of course, adopt extensively pa- 
rameterized models. But they also have the limitation of be- 
ing non-theoretical, i.e., not requiring physical insight they 
most often do not result in any. Indeed, simulations can be 
high-resolution experiments of undetermined accuracy for any 
physical system. 

Quasi-chemical theory is not a take-it-or-leave-it model, 
and not a series expansion, but a well-defined structure for 
combining computational results from distinct sources that 
treat separately near and more distant neighbors. As more ad- 
vanced applications are encountered — here with the Rb + (aq), 
which localizes near-neighbor water molecules slightly less 
definitely than some less advanced cases — some physical 
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learning and judgement is required. A big step in that learn- 
ing has been to catagorize near-neighbors on the basis of the 
neighborship decomposition of the radial distribution as in 
Fig. [2] This is in contrast to identification of neighbors on 
the basis of the location of the first minimum of that radial 
distribution, which is often less than compelling. Another 
step in that learning has been to consider more sophisticated 
procedures for estimation of cluster hydration free energies. 
Our discussion here has emphasized further the clear learning 
point that fluctuations of the composition of the inner-shell are 
numerically non-significant for strongly bound cases such as 
Rb + and where a well-informed identification of inner-shell 
ligands has been achieved. 
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